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ResearchObjectives 


The  overall  objectives  of  this  work  are  to: 

•  Quantify  the  requirements  for  controlled  precipitation  of  radiation  belt  particles  in  order  to  mitigate  space 
particle  effects 

•  Address  the  issue  of  coupling  of  space  based  antennas  to  the  surrounding  magnetoplasma  for  use  in  system 
design 

•  Determine  radiated  power  from  single  transmitting  element 

The  completion  of  these  goals  requires  the  design  of  antenna-in-plasma  codes  which  simulate  the  interaction 
of  radiating  antennas  and  the  magnetized  plasma  in  which  they  are  immersed.  Different  codes  have  been  de¬ 
veloped  that  solve  various  aspects  of  the  coupling  problem  and  that  are  tailored  to  the  particular  physics  in  the 
corresponding  region  of  the  coupling  environment. 


1  Introduction 

The  primary  focus  of  our  research  effort  under  AFOSR 
sponsorship  has  been  the  development  of  software  for 
the  purpose  of  studying  the  coupling  of  antennas  to  a 
magnetized  plasma  for  application  in  the  remediation  of 
energetic  electrons  from  the  Van  Allen  radiation  belts 
through  in-situ  VLF  wave  injection.  As  a  direct  result 
of  this  research,  we  have  developed  time  and  frequency 
domain  codes  for  this  purpose  in  order  to  gain  insight 
into  the  complex  nature  of  the  physical  processes  inher¬ 
ent  when  using  dipole  antennas  in  a  magnetized  plasma 
operating  at  VLF  frequencies. 

Electric  dipole  antennas  are  commonly  used  in 
space  plasmas  with  applications  that  range  from  radio 
frequency  probing  of  the  magnetosphere  to  plasma  di¬ 
agnostics  [Carpenter  and  Anderson,  1992,  Carpenter 


et  al.,  2003,  Bell  et  al.,  2004,  Platino  et  at.,  2005],  With 
the  recent  interest  in  the  use  of  these  antennas  for  the 
purpose  of  radiation  belt  remediation  (RBR)  [Inan  et  al., 
2003],  the  characterization  of  the  antenna-plasma  cou¬ 
pling  behavior  in  this  regime  is  of  primary  importance. 
Along  with  ground  based  transmitters  and  naturally  oc¬ 
curring  sources  of  VLF  radiation  such  as  lightning  and 
chorus  [Abel  and  Thorne,  1998],  it  has  been  shown  that 
in-situ  injection  of  VLF  electromagnetic  waves  can  sub¬ 
stantially  reduce  population  lifetimes  of  highly  energetic 
electrons  in  the  Earths  radiation  belts.  The  dominant 
mechanism  behind  the  precipitation  of  these  energetic 
particles  is  pitch  angle  diffusion  via  cyclotron  resonant 
wave-particle  interactions  [Albert,  2001]. 

Over  the  past  five  decades,  there  has  been  signif¬ 
icant  interest  in  the  study  of  antenna-plasma  interac¬ 
tions.  Some  of  the  first  work  in  the  field  was  performed 
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by  Balmain  [1964]  who  analytically  determined  the  input 
impedance  of  electrically  short  dipole  antennas  using  a 
quasi-electrostatic  approach  in  a  collisional  ionospheric 
plasma.  Subsequent  review  papers  by  Balmain  [1972, 
1979],  extended  this  analysis  to  a  variety  of  antennas 
operating  in  both  isotropic  and  anisotropic  plasmas  in¬ 
cluding  resonance  and  nonlinear  effects.  At  roughly  the 
same  time  Wang  and  Bell  [1969,  1970,  1972c,a,b]  de¬ 
rived  closed-form  expressions  for  the  terminal  proper¬ 
ties  and  radiation  characteristics  of  both  dipole  and  loop 
antennas  operating  at  ELF/VLF  using  a  fully  electromag¬ 
netic  approach.  Early  theoretical  work  often  made  sim¬ 
plifying  assumptions  in  order  to  provide  analytical  solu¬ 
tions  such  as  a  triangular  current  distribution  on  the  an¬ 
tenna  and  a  homogeneous  linearized  cold/warm  plasma 
environment. 


Whether  operating  as  transmitting  or  receiving  ele¬ 
ments,  electric  dipole  antennas  will  be  surrounded  by 
an  electrostatic  sheath  which  can  significantly  alter  the 
antenna  properties  (both  near  and  far  field)  relative  to 
its  freespace  equivalent.  For  receiving  purposes,  the 
sheath  will  be  on  the  order  of  a  few  Debye  lengths  and 
is  well  approximated  by  existing  analytical  theory.  How¬ 
ever,  when  used  for  applications  requiring  the  applica¬ 
tion  of  large  potentials  to  the  transmitting  element  far  in 
excess  of  the  surrounding  plasma  potential,  the  sheath 
is  highly  nonlinear  and  its  structure  is  unknown. 


The  emphasis  of  this  report  is  on  the  near  field  prop¬ 
erties  of  dipole  antennas  including  the  nonlinear  effects 
of  the  sheath.  This  paper  proceeds  with  a  theoreti¬ 
cal  background  outlining  the  various  models  used  for 
simulating  antenna-plasma  interactions.  We  begin  with 
the  most  simplified  approach  constituting  a  linearized 
plasma-fluid  development  coupled  with  Maxwell's  equa¬ 
tions  in  order  to  determine  the  impedance  and  cur¬ 
rent  distribution  on  dipole  antennas  in  the  absence  of 
the  sheath.  Then  we  increase  the  complexity  of  our 
model  by  including  finite  temperature  and  nonlinear  ef¬ 
fects  by  extending  the  number  of  fluid  moments  in  our 
approach  in  order  to  determine  the  effects  of  the  sheath 
on  the  near  field  properties  of  dipole  antennas.  The  ter¬ 
minal  properties  of  dipole  antennas  for  the  aforemen¬ 
tioned  cases  are  provided  and  a  discussion  of  the  rela¬ 
tive  importance  of  various  fluid  quantities  in  time-varying 
sheath  formation  is  given. 


2  Theoretical  Development 

Our  plasma  formulation  uses  a  macroscopic  multi-fluid 
approach  to  solve  for  the  coupling  between  dipole  an¬ 
tennas  and  the  plasma  in  which  they  operate.  The  fluid 
mode!  is  comprised  of  moments  of  the  Vlasov  equation 
for  each  particle  species  representing  the  time-evolution 
of  the  particle  distribution  function  /  =  /(r,  v,  t)  in  a  colli¬ 
sionless  plasma  where  r,  v  and  /  are  coordinates  in  con¬ 
figuration  space,  velocity  space  and  time  respectively. 
The  Vlasov  equation  given  by  Equation  1 . 


where  m  is  the  particle  mass,  and  F  =  ofE  +  vxB) 
is  the  Lorentz  force  with  charge  q ,  electric  field  E,  and 
magnetic  field  B.  Since  the  bulk  of  the  plasma  between 
L  =  2  and  L  =  3  is  virtually  cold  (~  2000°)  consist¬ 
ing  of  very  low  energy  particles  [Bezrukikh  et  al.,  2003] 
we  assume  our  initial  distribution  for  each  species  is  a 
Maxwellian  distribution  as  in  Calder  and  Laframboise 
[1990],  Calder  et  al.  [1993]  and  given  by  Equation  2  in  a 
three-dimensional  Cartesian  coordinate  system. 

(2> 

v  =  [(vx  -  Ul0){vy  -  Uy0)(v1!  -  1!.o)] 

where  N0  is  the  bulk  density,  v  =  vfvx.vy.v.,)  is  the  ve¬ 
locity  vector  represented  by  coordinates  v*,  vy  and  va 
in  velocity  space,  vxo,  Vyo  and  u20  are  the  drift  velocities 
in  each  respective  direction,  k  is  Boltzmann’s  constant, 
and  T  is  temperature.  In  our  simulations,  we  assume 
that  the  drift  velocities  ul0.  Uj,g  and  vtQ  are  initially  zero. 

2.1  Moments  of  Vlasov  Equation 

The  series  of  moments  that  comprise  our  fluid  model  are 
derived  by  integrating  the  Vlasov  equation  over  velocity 
space  in  powers  of  vn  where  n  is  the  order  of  the  mo¬ 
ment.  In  many  cases  including  ours,  it  is  convenient  to 
define  the  second  and  third  order  moments  which  spec¬ 
ify  the  stress  tensor  n  and  energy  flux  density  £  in  the 
rest  frame  of  the  species  using  the  relation  c  =  v  -  u 
where  c  is  the  random  velocity  due  to  the  thermal  mo¬ 
tions  of  the  particles,  u  is  the  flow  velocity,  and  v  defines 
the  velocity  distribution  of  the  particles  in  phase  space 
as  mentioned  earlier.  In  the  rest  frame,  n  and  £  be¬ 
come  the  pressure  P  and  heat  flux  density  Q  which  are 
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integrated  with  respect  to  the  Viasov  equation  in  pow¬ 
ers  of  [v  -  u]n.  The  quantities  comprising  the  first  four 
moments  given  by  Equations  3a-3d  which  are  found  in 
Chust  and  Belmont  [2006]  correspond  to  density,  mo¬ 
mentum,  pressure,  and  heat-flux  respectively  and  are 
repeated  here  for  convenience. 
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In  Equations  3c  and  3d,  nr  represents  the  gyrofre- 
quency  vector  along  the  magnetic  field,  R  is  the  4,h  order 
moment,  and  the  sym  symbolizes  a  symmetric  tensor.  It 
is  useful  to  note  that  with  each  moment,  the  rank  of  each 
tensor  increases  by  xn  where  x  is  the  number  of  dimen¬ 
sions  and  n  is  the  order  of  the  moment,  Therefore,  the 
density  equation  for  each  species  comprises  only  a  sin¬ 
gle  element  while  the  heat-flux  density  tensor  which  is 
of  rank  three  consists  of  27  individual  components.  This 
can  quickly  render  the  fluid  method  intractable  unless 
simplifying  assumptions  are  made  for  higher  moments. 

We  will  compare  three  different  fluid  closure  relations 
in  this  paper  for  simulating  the  coupling  of  the  antenna 
to  the  plasma.  The  first  is  the  cold  plasma  approxima¬ 
tion  in  which  we  assume  that  the  pressure  tensor  P  is 
zero  in  the  momentum  transport  equation.  We  also  as¬ 
sume  small  scale  perturbations  in  each  of  the  conserved 
quantities  such  that  Equations  3a-3b  can  be  linearized 
and  combined  into  a  single  equation  for  the  current  for 
each  species  given  by  Equation  4. 

~  +  vja  =  —  (qanaE  +  JQ  X  B0)  (4) 

where  E  is  the  electric  field,  J,  v ,  n,  q,  and  m  are 
the  current  density,  collision  frequency,  number  density, 
charge,  and  mass  of  species  a  and  B  is  the  background 
magnetic  field.  Equation  4  represents  a  simplified  ver¬ 
sion  of  the  generalized  Ohm's  Law. 

The  second  closure  relation  is  the  isothermal  approx¬ 


imation  which  provides  a  relationship  between  tempera¬ 
ture  and  pressure  given  by  P  =  nkT.  This  is  the  equa¬ 
tion  of  state  based  on  the  ideal  gas  law  for  thermody¬ 
namics  and  provides  a  mechanism  for  truncating  the  set 
of  infinite  moments  at  the  momentum  equation  assum¬ 
ing  a  diagonaf  pressure  tensor  and  for  the  case  in  which 
the  temperatures  parallel  and  perpendicular  to  the  mag¬ 
netic  field  are  the  same,  simplifies  to  an  isotropic  pres¬ 
sure.  The  third  closure  relation  utilizes  the  assumption 
of  adiabadicity  V  Q  =  0  which  essentially  means  that 
the  contribution  due  to  the  heat-flux  is  negligible.  How¬ 
ever,  in  this  case,  all  components  of  the  pressure  tensor 
are  kept  in  the  fluid  description.  This  assumption  as  with 
all  other  closure  relations  given  assumes  that  the  dis¬ 
tribution  function  remains  sufficiently  compact  about  the 
mean  (thermal  velocity)  so  that  the  first  few  moments 
provide  an  adequate  description  of  the  distribution.  In 
fact,  up  until  the  addition  of  the  heat-flux  moment  which 
allows  for  skew  about  the  mean  velocity  in  the  distribu¬ 
tion,  all  other  moments  assume  a  symmetry  about  the 
mean  velocity.  In  the  Results  sections,  we  will  com¬ 
pare  via  simulation  various  truncation  relations  for  both 
magnetized  and  unmagnetized  plasmas  with  respect  to 
sheath  formation,  whereas  Chust  and  Belmont  [2006] 
provide  a  more  analytical  justification  based  on  various 
time  and  spatial  scales  for  the  various  truncation  mecha¬ 
nisms  without  specifically  addressing  sheath  dynamics. 

In  our  initial  cold  plasma  simulation  work,  the  field 
quantities  are  solved  for  using  Maxwell’s  equations  given 
by  Equations  5a  -  5b. 


V  X  H  = 

N 

(5a) 

dH 

(5b) 

V  x  E  = 

where  E  and  H  are  the  wave  electric  and  magnetic 
fields  and  J  is  the  auxiliary  current  defined  in  Equation 
4.  As  we  will  show  in  the  Cold  Plasma  Results  section, 
the  current  distribution  is  triangular  for  all  cases  shown 
using  the  fully  electromagnetic  model.  With  the  inclu¬ 
sion  of  additional  moments  defining  the  warm  plasma 
model  in  order  to  study  sheath  dynamics,  we  utilize  Pois¬ 
son's  equation  in  order  to  ease  the  computational  bur¬ 
den.  Poisson's  equation  is  given  by: 
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V  •  E  =  — —  (6) 

where  e„  is  the  permittivity  of  free  space,  and  the  charge 
density  p  is  summed  across  all  species  «  where  a  -  2 
representing  only  electrons  and  protons  in  our  model. 
The  use  of  Poisson’s  equation  is  fully  justified  since  the 
sheath  is  primarily  electrostatic  in  character  and  the  use 
of  Poisson’s  equation  enforces  a  triangular  current  distri¬ 
bution  across  the  antenna  elements  which  we  will  show 
is  consistent  with  our  cold  plasma  results. 


3  Cold  Plasma  Results 

We  now  present  results  for  the  current  distributions  and 
input  impedance  of  electric  dipole  antennas  in  a  mag¬ 
netized  plasma.  Finite  Difference  Frequency  Domain 
(FDFD)  simulations  are  carried  out  for  dipole  antennas 
oriented  perpendicular  with  respect  to  the  ambient  static 
magnetic  field.  This  orientation  is  chosen  since  the  an¬ 
tenna  pattern  and  power  delivery  are  optimal  for  launch¬ 
ing  waves  parallel  to  the  static  magnetic  field  Wang  and 
Bell  [1972c, a].  Antennas  considered  for  the  purpose  of 
our  application  are  on  the  order  of  100  m  in  length  and 
up  to  20  cm  in  diameter. 

The  FDFD  method  is  well  suited  to  model  small  ge¬ 
ometries  with  respect  to  a  free  space  wavelength.  The 
antenna  itself  is  assumed  to  be  a  PEC  (Perfect  Electric 
Conductor)  and  the  current  distribution  along  the  length 
of  the  antenna  is  calculated  by  taking  a  line  integral  of 
the  frequency  domain  magnetic  field  components  encir¬ 
cling  each  wire  element  along  the  length  of  the  antenna. 
The  input  impedance  is  calculated  using: 


Vij)  (fE-dl)frtil 
1(f)  (§H-dl)ftt: 


where  the  field  quantities  are  already  in  the  frequency 
domain  per  use  of  the  FDFD  method.  (7)  represents  the 
ratio  of  the  complex  phasors  quantities  for  the  current 
and  voltage  at  the  terminals  of  the  antenna. 

For  the  purposes  of  our  simulation,  we  first  examine 
the  properties  of  a  100  m  long  electric  dipole  antenna 
in  a  cold  magnetized  plasma  operating  near  L  -  2  in 
the  magnetic  equatorial  plane.  We  consider  an  electron- 
proton  plasma  with  fpt.  =  586  kHz  and  fce  =  110  kFiz 
in  a  collisionless  environment.  The  z  directed  antenna 


is  located  in  the  center  of  the  space  and  is  20  cm  in 
diameter  which  corresponds  to  the  smallest  cell  size  in 
the  space.  The  magnetic  field  is  oriented  in  the  +y  di¬ 
rection  and  a  perfectly  matched  layer  (PML)  boundary 
condition  is  used  to  truncate  the  space  in  all  directions. 
The  dipole  antenna  is  excited  with  an  E~  hard  source  in 
the  gap  between  the  conducting  elements  with  a  value 
of  1  V/m  and  the  system  is  allowed  to  converge  with  a 
relative  residual  norm  of  10"6. 

One  of  the  primary  benefits  of  using  frequency  over 
time  domain  analysis  is  the  ability  to  use  a  different  mesh 
and  PML  configuration  for  each  simulation  run.  Though 
we  do  not  use  this  advantage  to  the  full  extent  available 
(a  different  configuration  for  each  frequency),  we  do  use 
a  different  mesh  and  PML  configuration  for  frequencies 
below  and  above  /lhr  for  which  the  propagation  char¬ 
acteristics  are  quite  different  as  shown  in  Figures  1  and 
2. 


Figure  1:  Refractive  index  surface  for  /u, n  <  f  < 


Figure  2:  Refractive  index  surface  for  /„  <  /  < 


In  Figures  1  and  2,  fc  represents  the  wave  number, 
vg  is  the  group  velocity  or  velocity  of  energy  flow  given 
by  the  normal  to  the  refractive  index  surface,  and  £fr,.s  is 
the  resonance  cone  angle  at  which  the  refractive  index 
tends  to  infinity.  The  refractive  index  surfaces  depicted  in 
Figures  1  and  2  are  functions  of  the  wave  normal  angle. 
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defined  as  the  angle  between  the  ^-vector  direction  and 
the  ambient  magnetic  field.  The  Gendrin  angle  is  the 
non-zero  wave  normal  angle  at  which  the  group  velocity 
is  parallel  with  the  static  magnetic  field  Gendrin  [1961]. 
The  Gendrin  angle  0g  is  illustrated  in  Figure  1  and  for 
high  ratios  of  is  given  approximately  by  the  relation 
Gendrin  [1961]"' 


COS  fig  RS  - - ,  (8) 

^ ce 

where  u  and  ujC£  are  the  angular  wave  frequency  and 
electron  gyrofrequencies  respectively. 

For  frequencies  /  >  /lhr.  there  exists  a  range  of 
^•-vectors  for  which  the  refractive  index  is  very  large  and 
tending  to  infinity  at  the  resonance  cone  angle  8res  as 
shown  in  Figure  1.  It  is  therefore  imperative  to  utilize 
much  smaller  cells  in  order  to  capture  these  tiny  wave¬ 
lengths  resulting  from  the  high  refractive  index  relative 
to  those  used  for  frequencies  /  <  /lhr-  For  frequencies 
below  /lhr-  the  resonance  cone  disappears  as  shown 
in  Figure  2  with  the  refractive  index  surface  being  closed 
and  possessing  a  maximum  of  n  «  600  at  directions  or¬ 
thogonal  to  the  background  magnetic  field.  The  refrac¬ 
tive  index  surface  for  /  <  /lhr  becomes  more  isotropic 
with  decreasing  frequency  and  thus  larger  cells  sizes 
may  be  used. 

As  stated  previously,  the  computational  space  is 
truncated  using  a  PML  absorbing  boundary  condition 
(ABC)  Berenger  [1994],  As  with  the  computational 
mesh,  the  PML  configuration  is  different  for  frequencies 
above  and  below  /lhr-  For  frequencies  /  >  /lhr.  the 
PML  consists  of  10  cells  in  the  x  direction,  15  cells  in 
the  i)  direction,  and  10  cells  in  the  z  direction.  The  PML 
layers  in  both  the  x  and  z  directions  are  made  to  ab¬ 
sorb  only  evanescent  waves,  while  the  PML  layers  in 
the  ij  direction  absorb  both  propagating  and  evanescent 
waves  in  this  frequency  range.  For  frequencies  below 
/lhr.  there  are  10  PML  layers  in  all  directions  and  each 
PML  is  made  to  absorb  both  propagating  and  evanes¬ 
cent  waves.  The  PML  parameters  are  different  for  fre¬ 
quencies  above  and  below  /lhr  since  the  cell  sizes  and 
refractive  index  surfaces  are  quite  different  in  the  two 
cases.  The  PML  performance  up  to  10  kHz  including 
frequencies  above  and  below  /lhr  for  the  simulations  in 
the  L  =  2  environment  are  shown  in  Figures  3  and  4 
corresponding  to  PML  orientations  parallel  and  perpen¬ 
dicular  to  the  static  magnetic  field  respectively. 


- RliCP  (Propagating)  -  -  -  LUCS*  iKvancsccnO  | 


Figure  3:  Reflection  coefficient  calculations  for  PML  oriented 
along  ./-direct ion,  parallel  to  the  static  magnetic  field,  for  an  an¬ 
gle  of  15°  with  respect  to  normal  incidence. 


"  RHCP  (Propagating) 

- RHCP  (Evanescent)  -  » -  LHCP  (Evanescent) 


Figure  4:  Reflection  coefficient  calculations  for  pml  oriented 
in  x  and  5 -directions,  perpendicular  to  the  static  magnetic  field, 
for  an  angle  of  45°  with  respect  to  normal  incidence. 


There  are  several  things  to  notice  about  the  plots  of 
Figures  3  and  4.  First  of  all,  the  only  propagating  modes 
in  the  frequency  range  /a  <£  /  <  fee  where  /„  <  /lhr 
are  right-hand  circularly-polarized  (RHCP).  All  waves 
launched  from  the  antenna  that  are  left-hand  circularly- 
polarized  (LHCP)  are  evanescent  in  the  plasma  at  these 
frequencies.  The  discontinuity  in  the  reflection  coeffi¬ 
cient  calculations  at  /lhr  is  a  direct  result  of  the  differ¬ 
ences  in  mesh  and  PML  geometries  across  this  transi¬ 
tion  region  as  stated  earlier.  Though  the  PML  perfor¬ 
mance  for  LHCP  waves  representing  the  directions  or¬ 
thogonal  to  the  static  magnetic  field  described  by  Figure 
4  is  relatively  poor,  these  waves  will  reflect  into  the  PML 
parallel  to  the  static  magnetic  field  of  Figure  3  and  be 
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absorbed  with  greater  attenuation.  The  PML  in  the  di¬ 
rection  orthogonal  to  the  static  magnetic  field  as  repre¬ 
sented  in  Figure  4  for  frequencies  /  >  /lhr  has  been 
tailored  to  absorb  evanescent  waves  only  to  avoid  the 
PML  instabilities  mentioned  earlier.  As  a  result,  the  in¬ 
cident  RHCP  propagating  modes  experience  no  attenu¬ 
ation  and  are  perfectly  reflected.  Finally,  for  small  inci¬ 
dent  angles,  the  RHCP  wave  incident  upon  the  PML  is 
evanescent.  The  resonance  cone  angle  of  Figure  t  is 
5°  at  9.2  kHz  and  marks  the  point  at  which  the  inci¬ 
dent  RHCP  waves  become  evanescent  and  able  to  be 
absorbed  by  this  PML. 


3.1  Current  Distributions  and  Input  Impedance 
Calculations  for  a  100  m  Antenna  at  L  =  2 


The  first  case  study  is  a  100  m  antenna  located  in  the 
equatorial  plane  at  L  =  2.  We  compare  the  current 
distributions  for  frequencies  above  and  below  the  local 
LHR  frequency  which  is  /l hr  -  2.55  kHz.  Figures  5  and 
6  represent  the  current  distributions  for  two  frequencies 
below  /LHr. 


Figure  5:  L  =  2  Current  distribution  (or  a  100  m  antenna  at  f  = 
400  Hz  driving  frequency. 


position  along  antenna  (m) 


Figure  6:  L  =  2  Current  distribution  (or  a  100  m  antenna  at  f  = 
2.0  kHz  driving  frequency. 

It  can  be  seen  from  Figures  5  and  6  that  the  current 
distributions  are  virtually  identical  to  the  assumed  trian¬ 
gular  distribution  of  Wang  and  Bell  [1969,  1970],  Wang 
[1970].  One  important  point  is  that  for  a  simulated  an¬ 
tenna  of  finite  thickness,  the  current  is  non-zero  at  the 
ends,  contrary  to  the  ideal  case,  since  the  finite  area 
allows  for  a  build  up  of  charge  at  the  tips.  The  simula¬ 
tion  results  thus  reflect  this  realistic  condition  much  bet¬ 
ter  than  the  idealized  case  shown  in  dashed  lines.  Fig¬ 
ures  7  and  8  represent  the  current  distributions  for  two 
frequencies  above  /lhr- 


Figure  7:  L  =  2  Current  distribution  for  a  100  m  antenna  at  f  = 
2.6  kHz. 
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Figure  8:  L  =  2  Current  distribution  for  a  100  m  antenna  at  f  = 
10,0  kHz. 

Once  again,  there  is  no  significant  deviation  from  the 
assumed  triangular  distribution,  except  for  the  realistic 
end-effect  due  to  the  finite  antenna  radius. 

Figure  9  compares  the  simulated  input  impedance  of 
the  100  m  dipole  antenna  at  L  =  2  with  results  obtained 
from  Wang  and  Bell  [1969, 1970],  Wang  [1970]. 


Figure  9:  Input  impedance  for  a  100  m  antenna  at  L  ~  2. 


The  plots  for  both  the  resistance  and  reactance  cal¬ 
culated  with  our  numerical  simulation  are  in  good  agree¬ 
ment  with  those  evaluated  analytically  by  Wang  and  Bell 
[1969,  1970],  Wang  [1970].  Below  /lhr.  Wang  and  Bell 
[1969,  1970],  Wang  [1970],  predict  the  reactance  to  vary 


from  approximately  lOOfi  at  zero  frequency  to  oo  at  the 
LHR  frequency.  Unlike  with  the  quasi-electrostatic  as¬ 
sumption  of  Balmain  [1964],  the  work  of  Wang  and  Bell 
[1969, 1970],  Wang  [1970]  predict  the  resistance  to  have 
a  non-zero  value  below  the  LHR  frequency  ranging  from 
OH  at  zero  frequency  to  oo  at  the  LHR  frequency.  These 
trends  are  reflected  in  the  simulated  results  as  shown  in 
Figure  9.  Above  /LHr.  the  analytical  reactance  varies 
from  Oil  to  about  -kil  at  10  kHz.  The  simulated  results 
in  these  regimes  are  within  about  15fi. 


3.2  Current  Distributions  and  input  Impedance 
Calculations  for  a  100  m  Antenna  at  L  =  3 

The  second  cold  plasma  case  study  examines  the  prop¬ 
erties  of  a  1 00  m  antenna  located  at  L  =  3  in  the  equato¬ 
rial  plane.  Typical  values  of  the  plasma  and  gyrofrequen- 
cies  at  L  =  3  are  /pG  =  220  kHz  and  fce  =  33  kHz  respec¬ 
tively.  Since  the  PML  performance  characteristics  of  Fig¬ 
ures  3  and  4  are  very  similar  to  the  simulation  setup  at 
L  =  2,  these  characteristics  are  not  shown.  Oniy  two 
different  examples  of  the  current  distribution  are  given 
since  they  are  not  markedly  different  than  those  for  the 
cases  of  an  antenna  at  L  =  2.  Figures  10  and  1 1  repre¬ 
sent  the  current  distributions  for  frequencies  above  and 
below  /uhr  —  761  Hz. 


position  along  antenna  (m) 


Figure  10:  L  =  3  Current  distribution  for  a  100  m  antenna  at 
f=460  Hz  driving  frequency. 
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Figure  11:  L  =  :i  Current  distribution  for  a  100  m  antenna  at 
f =4.0  kHz  driving  frequency. 

It  can  be  seen  from  Figures  10  and  1 1  that  there  is 
no  significant  deviation  from  a  triangular  current  distri¬ 
bution  for  an  antenna  subject  to  a  decrease  in  plasma 
and  gyrofrequencies.  Figure  12  compares  the  simulated 
input  impedance  of  the  100  m  dipole  antenna  at  L  =  3 
with  results  obtained  from  Wang  and  Bell  [1969,  1970], 
Wang  [1970], 


Figure  12:  Input  impedance  for  a  100  m  antenna  at  L  =  3. 


As  with  the  results  from  the  previous  case  study,  the 
impedance  characteristics  in  Figure  12  exhibit  very  good 
agreement  with  the  work  of  Wang  and  Bell  [1 969, 1 970], 
Wang  [1970].  Our  results  from  this  section  validate  our 


model  using  the  cold  plasma  approximation.  In  the  next 
section,  we  extend  the  capability  of  the  model  through 
the  addition  of  more  moments  to  the  fluid  representation 
allowing  for  the  formation  of  the  sheath  surrounding  the 
antenna.  In  particular,  we  will  examine  the  effects  of  the 
sheath  on  the  impedance  characteristics  of  the  antenna. 


4  Warm  Plasma  Results 

The  warm  plasma  fluid  formulation  allows  for  the  study 
of  the  sheath  formation  surrounding  an  electric  dipole 
antenna.  We  will  show  that  it  is  the  sheath  which  domi¬ 
nates  the  tuning  properties  of  the  antenna  and  it  is  thus 
imperative  to  accurately  simulate  its  contribution  to  the 
overall  antenna-plasma  coupling  response.  In  this  sec¬ 
tion,  we  will  examine  the  effects  of  two  various  closure 
mechanisms  in  the  context  of  2-dimensional  simulations 
and  then  analyze  the  sheath  surrounding  an  electric 
dipole  antenna  in  3-dimensions, 

4.1  Closure  Relations 

In  the  warm  fluid  approximation,  no  assumptions  of  lin¬ 
earity  are  made  and  thus  all  convective  terms  in  the  sys¬ 
tem  of  moments  defined  by  Equations  3a-3d  are  pre¬ 
served,  Two  warm  plasma  closure  mechanisms  are  con¬ 
sidered  in  this  report.  The  first  closure  relation  is  the 
isothermal  approximation  which  assumes  that  the  pres¬ 
sure  tensor  is  diagonal  and  given  by  the  ideal  gas  law 
relation  discussed  earlier.  In  this  approximation,  the  sys¬ 
tem  of  equations  representing  the  plasma  is  given  by  the 
first  two  fluid  moments  given  by  Equations  3a-3b. 

The  second  closure  relation  assumes  that  the  sys¬ 
tem  is  adiabatic  V-Q  =  0  and  thus  all  terms  involving  the 
heat  flux  tensor  are  zero.  In  this  formulation,  the  plasma 
is  represented  by  the  first  three  moments  of  Equation  1 
given  by  Equations  3a-3c. 

4.2  Two-dimensional  Simulation 

The  first  study  considered  is  that  of  a  two  dimensional  in¬ 
finite  line  source  in  an  isotropic  plasma  with  parameters 
corresponding  to  L  =  2.  We  use  an  artificial  ion-electron 
mass  ratio  of  ^  =  200  which  is  a  commonly  used  prac¬ 
tice  in  order  to  ease  the  computation  burden  [Calder 
et  al.,  1993].  The  simulation  uses  Neumann  boundary 
conditions  for  the  fluid  at  the  edge  of  the  space  and  any 
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particles  hitting  the  line  source  are  completely  absorbed. 
A  plot  of  the  simulation  domain  is  shown  in  Figure  13. 


Figure  13:  2-dimensional  computational  domain. 


We  examine  the  time-varying  sheath  formation  us¬ 
ing  a  quasi-electrostatic  approach  as  discussed  earlier. 
This  involves  the  use  of  Poisson’s  equation  in  which  a 
sinusoidal  potential  is  applied  to  the  line  element  with  a 
peak  potential  of  approximately  1 00V.  In  order  to  deter¬ 
mine  the  importance  of  the  proton  dynamics  in  sheath 
formation,  we  consider  RF  frequencies  above  and  be¬ 
low  the  proton  plasma  frequency.  Since  the  plasma  is 
isotropic  in  our  2-d  case,  we  take  a  vertical  slice  plane 
through  the  line  element  and  plot  the  resulting  potential 
and  density  profiles  using  a  1-d  representation.  A  dia¬ 
gram  of  the  slice  plane  is  shown  in  Figure  14. 


BJ 


Figure  14:  Slice  plane  through  infinite  line  source. 


As  the  potential  becomes  positive,  the  protons  are 
pushed  away  from  the  antenna  white  the  electrons  are 
collected,  while  during  the  negative  cycle,  the  opposite 
occurs  with  the  protons  being  collected  and  the  elec¬ 
trons  being  repelled  from  the  line  element.  Figures  15 


and  16  show  the  potential  and  density  profiles  for  the 
electrons  and  protons  of  our  two  species  plasma  after  a 
number  of  periods  of  the  25kHz  applied  potential  {5kHz 
above  the  proton  plasma  frequency)  by  which  time  the 
response  is  periodic. 
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Figure  15:  Snapshot  during  negative  potential  cycle  of  infinite 
line  source. 
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Figure  16:  Snapshot  during  positive  potential  cycle  of  Infinite 
line  source. 


During  the  start  of  the  initial  transient  response,  the 
protons  are  pushed  away  from  the  line  source  and  due 
to  their  inertia,  are  never  pulled  back  during  the  oppo¬ 
site  cycle.  This  forms  a  permanent  evacuated  region 
of  protons  with  the  only  substantial  current  collection 
on  the  line  element  being  due  to  the  electron  response. 
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This  behavior  is  in  disagreement  with  previous  analytical 
works  in  which  the  proton  motion  was  neglected  since  it 
was  thought  that  the  much  heaver  mass  of  the  protons 
would  essentially  make  them  immobile  [Mlodnosky  and 
Garriott,  1963,  Baker  et  al.,  1973,  Song  et  al.,  2007], 
Figures  1 7-20  represent  the  current  collection  on  the 
line  element  for  a  given  RF  applied  potential  for  fre¬ 
quencies  above  and  below  the  proton  plasma  frequency. 
The  blue  curves  represent  the  isothermal  approxima¬ 
tion  (2-moments)  and  the  red  curves  represent  the  adi¬ 
abatic  approximation  (3-moments).  The  proton  plasma 
frequency  with  the  proton-electron  mass  ratio  of  200  is 
approximately  20kHz. 
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Figure  17:  Current-voltage  relationship  for  isotropic  plasma 
with  sinusoidal  applied  potential  at  15kHz.  'Blue'  -  2  moments. 
Red'  -  3  moments. 
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Figure  18:  Current-voltage  relationship  for  isotropic  plasma 
with  sinusoidal  applied  potential  at  25kHz.  'Blue1  -  2  moments. 
’Red’  -  3  moments. 


Figures  17  and  18  represent  the  current  collection 
on  the  line  element  assuming  an  isotropic  plasma.  The 


current  waveform  is  shown  in  both  linear  and  log  scales 
in  order  to  visualize  the  small  contribution  of  the  pro¬ 
tons.  Although  the  proton  contribution  to  the  current  is 
negligible  as  seen  per  the  diode  type  behavior,  it  is  non¬ 
zero.  The  other  important  point  is  that  there  is  little  differ¬ 
ence  between  the  isothermal  and  adiabatic  approxima¬ 
tions  suggesting  that  only  a  small  number  of  moments 
are  necessary  in  order  to  capture  the  relevant  physics  in 
sheath  formation  even  under  the  collisionless  assump¬ 
tion.  This  conclusion  is  supported  by  the  earlier  work  of 
Thiemann  et  a!.  [1992]  who  compared  both  a  2-moment 
fluid  code  and  a  particle  in  ceil  (PIC)  code  with  good 
qualitative  agreement  in  a  collisionless  isotropic  plasma. 
Figures  19  and  20  include  the  effects  of  a  static  mag¬ 
netic  field  oriented  in  the  i  direction  as  shown  in  Figure 
13. 


Figure  19:  Current-voltage  relationship  for  magnetized  plasma 
with  sinusoida!  applied  potential  at  15kHz.  Blue'  -  2  moments. 
Red'  -  3  moments. 
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Figure  20:  Current-voltage  relationship  for  magnetized  plasma 
with  sinusoidal  applied  potential  at  25kHz.  'Blue'  -  2  moments. 
'Red'  -  3  moments. 
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The  inclusion  of  the  magnetic  field  in  the  formulation 
has  the  effect  of  limiting  the  current  flow  to  the  line  ele¬ 
ment  as  expected  since  charged  particles  gyrate  around 
magnetic  field  lines. 


In  this  paper,  we  define  the  edge  of  the  plasma 
sheath  as  the  point  at  which  the  potential  drops  to  within 
10%  of  the  maximum/minimum  amplitude  as  shown  in 
Figure  22. 


4.3  Three-dimensional  Simulation 

In  this  section  we  use  our  warm  fluid  code  to  examine 
the  formation  of  the  sheath  surrounding  a  dipole  antenna 
in  a  magnetized  plasma  in  3-dimensions.  The  dipole  is 
20  m  in  length,  1 0  cm  in  diameter,  and  possesses  a  2  m 
gap  between  elements.  We  consider  an  electric  dipole 
antenna  operating  in  conditions  corresponding  to  L  =  2 
and  L  =  a  in  the  geomagnetic  equatorial  plane  excited 
with  a  maximum  potential  difference  of  approximately 
100V.  The  same  scaled  proton-electron  mass  ratio  is 
used  as  in  the  2-dimensional  simulations.  Only  perpen¬ 
dicular  orientations  with  respect  to  the  background  mag¬ 
netic  field  are  considered  as  in  the  cold  plasma  simula¬ 
tion  runs.  Since  there  was  virtually  no  difference  in  the 
behavior  of  the  current  collection  for  the  2-dimensional 
cases  for  frequencies  just  above  and  below  the  proton 
plasma  frequency,  all  3-d  simulation  runs  are  performed 
at  the  proton  plasma  frequency.  We  assume  that  an 
electron  gun  is  present  on  board  the  spacecraft  such 
that  any  particles  hitting  the  antenna  are  completely  ab¬ 
sorbed  and  are  assumed  not  to  reverse  bias  the  an¬ 
tenna. 

The  first  case  study  is  a  dipole  antenna  at  L  =  3 
where  the  magnetic  field  is  not  as  strong  and  the  De¬ 
bye  length  is  slightly  smaller  than  the  operating  region 
at  L  =  2.  A  sketch  of  the  simulation  region  is  shown  in 
Figure  21. 

a  B0x 
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Figure  22:  Diagram  describing  location  of  sheath  edge. 


Figure  23  represents  and  orthographic  projection  of 
the  fully  expanded  sheath  during  the  peak  of  the  sinu¬ 
soidal  cycle  after  several  periods  in  each  slice  plane.  It 
can  be  seen  from  Figure  23  that  the  sheath  is  approxi¬ 
mately  t  m  wide  at  is  virtually  symmetric  in  each  of  the 
slice  planes  suggesting  that  the  magnetic  field  does  not 
play  much  of  a  role  in  the  steady-state  formation  of  the 
sheath.  In  addition,  the  sheaths  on  each  element  do  not 
overlap  meaning  that  both  elements  are  electrically  iso¬ 
lated  from  one  another  with  no  current  flowing  between 
them. 


Figure  21:  3-dimensional  computational  domain  at  L  =  3. 
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Figure  23:  Orthographic  projection  of  expanded  sheath  for  L  =  3  at  maximum  potential  difference  100V.  Operating  frequency  is 
at  /  = 


Figure  24  shows  the  time-domain  impedance  char¬ 
acteristics  for  the  dipoie  antenna  located  at  L  =  3  in¬ 
cluding,  potential,  current,  resistance  and  capacitance. 
The  most  notable  feature  of  Figure  24  is  that  both  the  re¬ 
sistance  and  capacitance  do  not  vary  significantly  over 
the  RF  cycle  with  the  capacitance  and  resistance  being 
roughly  10  pF  and  Ififcfi  respectively. 


Time  4is) 


Figure  24:  Time-domain  terminal  characteristics  for  20  m  an¬ 
tenna  located  at  L  =  3. 

The  final  case  study  is  a  dipole  antenna  located  at 
L  -  2  where  the  magnetic  field  is  more  appreciable  than 
that  which  exists  at  L  =  3.  Since  the  Debye  length  is 
smaller  at  L  =  2,  the  computational  domain  is  slightly 
smaller.  A  sketch  of  the  simulation  region  is  shown  in 
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Figure  25. 


BJ 


Figure  25:  3-dimensional  computational  domain  at  L  =  2. 


Although  not  shown  here,  the  only  major  difference 
between  the  runs  at  L  =  2  and  L  =  3  is  during  the 
initial  transient  response  in  which  the  transient  sheath 
expands  along  the  field  line.  However  the  steady-state 
sheath  (reached  after  a  few  sinusoid  cycles)  is  virtually 
the  same  as  that  found  at  L  —  3  with  only  a  very  slight 
increase  in  the  amount  of  asymmetry  along  the  slice 
planes.  Figure  26  represents  the  corresponding  ortho¬ 
graphic  projection  of  the  fully  expanded  sheath  for  the 
case  at  L  =  2. 


Figure  26:  Orthographic  projection  of  expanded  sheath  for  L  =  2  at  maximum  potential  difference  100V.  Operating  frequency  is 
at  /  = 


Figure  27  is  a  plot  of  the  terminal  properties  for  the 
antenna  at  L  =  2  with  similar  agreement  to  the  case 
study  at  L  =  3. 
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Figure  27:  Time-domain  terminal  characteristics  for  20  m  an¬ 
tenna  located  at  L  =  2. 


5  Discussion 

This  work  constitutes  a  major  achievement  in  the  de¬ 
termination  of  the  coupling  behavior  of  antennas  to  a 
magnetoplasma  using  a  multi  step  approach.  The  cold 
plasma  model  was  the  simplest  approximation  used  and 
from  the  results  we  were  able  to  obtain  some  insight 
into  the  coupling  of  dipole  antennas  to  the  surrounding 
plasma.  From  the  cold  plasma  simulation  results,  we 
determined  that  the  current  distribution  was  triangular 
for  the  antenna  parameters  used  which  validates  the  as¬ 
sumption  used  in  previous  analytical  work.  Also,  in  the 
absence  of  the  sheath,  it  was  shown  that  for  frequencies 
above  the  LHR  frequency,  the  input  impedance  does 
not  vary  significantly  as  a  function  of  frequency  mean¬ 
ing  that  the  antenna  is  self  tuning  and  can  be  used  over 
a  broad  frequency  range.  With  the  addition  of  the  warm 
plasma  fluid  approximation,  we  gained  the  capability  of 
examining  the  formation  of  the  electrostatic  sheath  and 
were  able  to  draw  a  number  of  conclusions  based  on 
these  results.  The  first  result  is  that  sheath  structure  is 
periodic  using  a  sinusoidal  waveform  excitation  exhibit¬ 
ing  a  quasi-steady  state  structure.  In  addition,  we  have 
shown  that  the  common  assumption  of  immobile  protons 
used  in  past  work  is  incorrect  with  the  density  of  pro¬ 
tons  varying  significantly  throughout  the  sheath  region 
and  contributing  slightly  to  the  current  collection.  Lastly, 
we  have  shown  that  not  only  does  the  sheath  dominate 
the  tuning  properties  of  the  antenna,  but  that  the  time- 
varying  resistance  and  capacitance  throughout  the  RF 


cycle  do  not  vary  by  orders  of  magnitude  (logarithmi¬ 
cally)  as  suggested  by  previous  authors  [Mlodnosky  and 
Garriott,  1963,  Baker  et  al.,  1973,  Song  et  al.,  2007], 
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